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Abstract 

A massively recurrent neural network responds on one side to input 
stimuli and is autonomously active, on the other side, in the absence of 
sensory inputs. Stimuli and information processing depends crucially 
on the qualia of the autonomous-state dynamics of the ongoing neural 
activity. This default neural activity may be dynamically structured 
in time and space, showing regular, synchronized, bursting or chaotic 
activity patterns. 

We study the influence of non-synaptic plasticity on the default dy- 
namical state of recurrent neural networks. The non-synaptic adaption 
considered acts on intrinsic neural parameters, such as the threshold 
and the gain, and is driven by the optimization of the information en- 
tropy. We observe, in the presence of the intrinsic adaptation processes, 
three distinct and globally attracting dynamical regimes, a regular syn- 
chronized, an overall chaotic and an intermittent bursting regime. The 



intermittent bursting regime is characterized by intervals of regular 
flows, which are quite insensitive to external stimuli, interseeded by 
chaotic bursts which respond sensitively to input signals. We discuss 
these finding in the context of self-organized information processing 
and critical brain dynamics. 



1 Introduction 



In the last couple of decades self organized processes have attracted the interests 
of many researchers from various scientific areas, both in natural and in social 
sciences. A system is said to be self-organizing, quite generally, when a state 
of high dyna mical complex i ty arises reliably from relatively simple basic orga- 



nizat i on ru l es (|Ashbv 



2003 



Grosl . l2010l ). 



1962 



Camazine. Deneubourg. Franks. Sneyd. fc Theraulal . 



It is often the case that the self-organization in dynamical systems is achieved 
through an interplay or regulative forces involving positive and negative feed- 
back, viz. through the interplay of internal drives which act destabilizing and 
regulating, respectively, onto the dynamics of the system. In general one type 
of feedback can dominate, driving the system towards a chaotic or towards an 
ordered phase respectively. A proper balance of the two opposing drives can bring 
the dynamical state at the point of a phase transition, a critical state. One speaks 
of self-organized criticality (SOC) whenever this balance is not achieved through 



the actions of an outside contro 



(IBak. Tang, fc Wiesenfeld 



1987 



l er bu t through internal self-organizing proces ses 



1998 



Bak. fc Paczuski 



1995 



Adami 



19951) 



As a dynamical system approaches a critical point, its spatiotemporal com- 
plexity rises. It has been suggested that this rise in the complexity improves 
the computational properties and the c apabi l ity of dynamical systems to process 



information (e.g.. ISqle. fc Miramontes 



Legenstein. fc Maass 



1995 



Bertschinger. fc Natschlager 



2004 



20061 ). This notion of "Computation at the edge of chaos' 



may als o been seen i n the broader context of "Life at the edge of chaos" (see 



Zimmer 



1999 



Gros . 



2010l ): the dynamical systems at the underpinning of all 



living have the tendency to self-organize themselves close to a critical state. 
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In recent years there have been many studies of the possible occurrence of SOC 
in neural networks with synaptic plasticity. Most of these studies have concluded 
that synaptic plasticity drives the dynarnics generically far below a critical point 



Siri. Quov. Delord. Cessac. fc Berrv 



2007 



Siri. Berry. Cessac. Delord. fc Quoyl . 



2008; 



Dauce. Quoy. Cessac. Doyon. &: Samuelided . Il998l ) , viz. it over- regulates the 



network dynamics. Hence, an organizational principle is needed which will main- 
tain an intermediate level of excitability in neural networks, preventing the occur- 
rence of dynamical states which are non- or hyper-reactive to external influences. 

It has been assumed for many years that the dominant driving force, shaping 
the brain's dynamical state, is the synaptic plasticity. Thus, little attention was 
put to ot her forms of neural adaption , the non-synaptic adaptation of individual 
neurons (IMozzachiodi. fc Byrnd . l2010l ). also known as intrinsic plasticity. Intrinsic 
plasticity is mostly manifested as a change in the excitability of a neuron, where 
this change is achieved through the adaptation on the level of membrane compo- 
nents. Here, we investigate the role of non-synaptic plasticity in the formation of 
complex patterns of neural activity. 



W e study a previously propo sed model of intrinsic plasticity (see 



Triesch 



2005 



2007 



Stemmler. fc Koch 



19991 ). and its influence on the dynamical properties of 



autonomous recurrent neural networks with rate encoding neurons in discrete time. 
Within this model neurons aspire to achieve, as an average over ti me, a firing - 
rate distribution which maximizes the Shannon information entropy ( iGrosl . l2010l ) . 
Therefore, the neurons are trying to homeostatically regulate a n entire distribution 
funct ion, a mechanism denoted polyhomeostatic optimization ( iMarkovic. fc Grod . 
20101 ). 

Intrinsic plasticity in the form of polyhomeostatic optimization gives rise, for 
random recurrent network topologies, to ongoing and self-sustained neural activi- 
ties with non-trivial dynamical states. Depending on the target mean firing rate 
and network parameters, one observes three distinct phase states: synchronized 
oscillations, an intermittent-bursting and a chaotic phase; all states being globally 
attracting in their respective phase spaces. 

In Section [2] we derived the stochastic learning rules for intrinsic adaptation. 
This is followed by the analysis of a single self-coupled neuron with intrinsic plas- 
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ticity (Section [3]) and the analysis of recurrent network with intrinsic plasticity 
(Section H]). Concluding remarks and discussion are finally provided in Section [H 



2 Stochastic adaptation 

We used a basic discrete-time, rate-encoding artificial neuron model. The firing 
rate y G [0, 1] of the neuron is given as a nonlinear transformation of the total 
synaptic input current x G (— oo,+oo), y = g{x). The transfer function g has a 
sigmoidal form, a usual choice being the logistic function 

y{t + 1) = ga,b{x{t)), ga,b{z) = ^ ^ l_az-b ' (1) 

where a is the gain and h the bias. The parameters of the transfer function, intrinsic 
parameters, will eventually become slow variables with a stochastic learning rules 
determining their time evolution. 

Let us denote with Px{x) the probability density function (PDF) of the total 
input. Given the relation ^ between the input current x and the output activity 
y we find 

Pa,b{y) = ^{y- 9a,bix)) Px{x) dx = ^ \x=g-\{y) (2) 

J-oo 9a,bW 

for Pa,b{y), the PDF of the firing rate. The main idea behind the derivation of 
adaption rules for the intrinsic parameters a and b is the assumption that the 
neuron's excitability should change in a way which maximizes the entropy of the 
firing rate dis tribution Vn^bi y), keeping at the same time the average activity at a 
desired level (ITrieschl . l2005l ) . 

The rational for this procedure is the following: the maximization of the firing 
rate entropy implies that a neuron will use the entire range of available activity 
states, optimizing the information transfer between neural input and output. Fur- 
thermore, the regulation of the average firing rate is present due to environmental 
constraints on the neuron, e.g. the limited energy resources needed for metabolic 
processes. 

Having a positive-definite variable y, with a fixed first moment, the maximum 



4 



entropy PDF corresponds to the exponential distribution 

1 



pxiy) 



Z{X) 



-Xy 



y e [0, 1] 



(3) 



where Z{\) = e ^^dy is the partition function. We will refer to the first moment 
of the PDF ([3]), denoted as /i, as the target average firing rate. It is given by 



A 



(4) 



In general the inverse function A(/i) cannot be found, but for A ^ 1 we recover 
the A — 7- which is valid for exponential PDFs defined on [0, oo]. 

A natural way to introduce a distance m easure between two PDFs is the 
Kullback-Leibler (KL) divergence (iGrod . l2010l ). defined as 



Dx{a,b) 



Pa,b{y) In ( ^"'f^.M dy 
V Pxiy) J 

-Er>A^^g'a,,{x)] + \E,^[gaA^)] - H\p,] + InZ(A) 



(5) 



where Ep^ [■] denotes the expectation value with respect to Px{,x), and H[p^ denotes 
the differential entropy functional. By minimizing Dx{a,b) with respect to the 
intrinsic parameters a and b one obtains the learning rules. In Eq. only the 
first two terms are functions of a and b. The gradient descent hence gives the 
following relation 



dDx{a,b) 
da 



E„ 



^ln<,(x)-AA,„,(x) 



Px{x) 



dx , 



(6) 



with a E {a,b}. As the input distribution Px{x) is in general unknown , it is 
suitable to derive the adaption rules by using a stochastic gradient descent (ISpalll . 
20051 ). Such adaption rules, for the update of the internal parameters a and b, are 
obtained by using the expression between the brackets on the right-hand side of 
Eq. (E]) . An advantageous side effect of this approach is that the adaptation rules 
become local in time. Using Eq. ([1]) for the transfer function ga,b{x) to evaluate 
Aa and A^, we obtain the stochastic learning rules 
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a{t + l) = a{t)+eAa{t) = a{t) + e {l/a{t) + x{t)A{t)) 

6(t + l) = bit)+eAt{t) = 6(t)+eA(t), (7) 



where e is the learning rate and 



A(t) 



{2 + \)yit + l) + \yit + lY . 



The learning rate e is assumed to be small; viz the time evolution of the internal 
parameters is slow compared to the evolution of both x{t) and y{t). In this way the 
stochastic adaptation, which depends only on the instantaneous values of the vari- 
ables, can closely match the direction of the deterministic gradient Eq. (|6]). Also, 
the input distribution can, in general, be non- stationary, therefore the minimum 
of the cost function Dx{a, b) could vary in time. For this reason, the learning rate 
should also be large enough for the adaptation to follow the changing minimum. 
However, any finite and constant learning rate e > doesn't satisfy the condition 
for exact co nvergence o f internal parameters into the minimum of Kullback-Leibler 
divergence (ISpalll . l2005l ). Still, if the learning rate decreases with every time step, a 
condition needed for strict convergence, the intrinsic adaptation will react slower 
with time to a variability in the position of the minimum of the cost function. 
Thus, it's more favorable to have here a constant learning rate which will result in 
the oscillations of the intrinsic parame ters around th e minimum. The amplitude 
of this oscillations roughly scales as e (iBottoul . 120041 ). thus a small learning rate 
also ensures convergence within a small vicinity from the minimum of the cost 
function D\{a, h). 



3 Single neuron 

We analyze initially a minimal network setup, a self-coupled neuron adapting 
homeostatically the intrinsic parameters of the transfer function. A synaptic con- 
nection between the axon and the dendrites of the same neuron is also known as 
an autapse. 
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Figure 1: (Left) A dependence of relative parameter change A (Eq. [9]) on output 
activity y for different target firing rates /x. (Right) Critical value of the gain Oc, 
as a function of the average firing rate /i. The colored area shows the region of 
stability of the fixpoint (?/*(A), 6*(a, A)), see Section[3l 

Neurons with an autapse are not rare in the brain. They have been observed 
in various brain regions and in different types of neurons. The d iscovery o f func - 
tional autapses pro v ides c lues for possible physiological roles ( jBekkersl . |2003| ). 
Herrmann, fc Klaud (120041 ) suggest that autapses lead to oscillatory behavior in 
otherwise non-oscillating neurons. We shall see below how this type of behav- 
ior spontaneously arises in self-excitatory neurons with intrinsic plasticity. We 
focused on the analysis of excitatory autapse and used it as a basis for under- 



standing the observed behavior in a larger network set up (see Section 



1). Fo r the 



20101). 



case of self-inhibition please refer to a separate study ( iMarkovic et al. 

The autapse neuron is equivalent to the identification x — ?/ in Eqs. ([1]) and 
(171). The complete set of evolution rules for the dynamical variables y{t), a{t) and 
h{t) is then 



y(^ + l) = ga(t)Mt)iyit)) 

bit + 1) = bit) + eA{t) 

a(t + l) = a{t) + e{l/a{t) +y{t)A{t)) 



(8) 



with 



A(t) = 1 - (2 + X)yit + 1) + \y^{t + 1) . (9) 

The right-hand side of ([9]) depends directly on y{t + 1) and only implicitly on y{t), 
as one can easily verify when going through the derivation of the rules ([7]) for the 



v[t) 

ait) 

m 




Figure 2: (Top) The time dependence of firing-rate y{t) (solid line), bias h{t) 
(squares) and gain a{t) (circles) for the one-site problem ([H]), with a learning rate 
e = 0.01 and a target average firing rate /i = 0.25; The gain a{t) is set initially 
bellow a critical value and since e ^ 1 the system relaxes quickly to the fixpoint of 
y = ga,b{y)- Once a{t) surpasses a certain threshold, compare Fig. [H the fixpoint 
becomes unstable and the system follows a limiting cycle. 

(Bottom) Maximal local Lyapunov Xmax{t) exponent compared to a Lyapunov 
exponent of a perturbation parallel to the flow A||(t). They were estimated along 
the points of the trajectory {y{t),a(t),b(t)). 

intrinsic plasticity. In plot at the left side of Fig. [1] we showed A(y) (Eq. ([9])) 
for various target firing rates fi. Note that A{y = 0) = 1 and A{y = 1) = — 1, 
independently of A. 



3.1 Stability analysis 

We first analyze a reduced model of the three evolution equations (|8]), obtained 
by setting Aa{t) = 0, viz considering a constant a = a(t) = a(t + 1). The reduced 
system contains a fixpoint (i/*(A), 6*(A, a)), where y* = [{2 + X) - a/4 + X'^]/2X and 
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b* = \n{y*/ (1 — y*)) — ay*. This fixpoint defines a one dimensional manifold in the 
complete phase space {y, b, a), where the stability of the manifold depends directly 
on the given value of a and A (see Fig. [1]). For a < Qc the dynamics is attracted 
toward the fixpoint {y*,b*), while for a > the fixpoint becomes repelling and 
the activity of the neuron follows a limiting cycle. One can show that the critical 
gain is = - y*). 

The time evolution of the full set of equations (see Fig. |2l top) approaches 
a limiting cycle, for all starting values of {y,b,a). The evolution rules ([8]) have 
fixpoint solutions also for a vanishing adapt ation, otz. fo r e — >^ 0. These fixpoints 



are turned for e > into attractor relics fiGrosl . 120071 . |2009| ). The trajectory 
slows down close to the attractor relics, giving rise to the transient firing states, 
observable in Fig. [21 This non-trivial activity pattern is a direct consequence 
of the polyhomeostatic adaption principle. The system cannot achieve, as an 
average over time, a non-trivial firing-rate distribution by settling into a steady 
state. Polyhomeostatic adaption hence forces the neuron to remain autonomously 
active, with varying firing rates. 

For an insight on the influence of external input signal on the dynamics, we have 
estimated the maximal local Lyapunov exponent Xmax{t) and the Lyapunov expo- 
nent for a perturbation in the direction of the flow A|| {Sf{t) oc {Ay{t), Aa{t), Ab(t))). 
They are presented in the bottom graph of Fig. [2l We see that the neuron is most 
sensitive to a perturbation during the transition between two attractor relics (low 
and high activity levels), since Xmax is positive through these transition periods. 
Also, A|| ~ Xmax during the fast transition between the attractor relics. We thus 
conclude that the direction of maximal sensitivity to perturbations is aligned with 
the direction of the flow at this points. Note that the two attractor relics can be 
stable during the same time period, although the activity settles in only one of 
them. This means that a transition could be induced with a sufficiently strong 
perturbation. 
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Figure 3: Output distributions of the two neurons with highest (diamonds) and 
lowest (circles) Kullback-Leibler divergence compared to the mean output dis- 
tribution (dashed line) and the target exponential output distribution (full line). 
The network size = 500 neurons and a target mean firing rate /i = 0.28 (top) 
and fi = 0.5 (bottom). The values are identical for all the neurons in the network 
and fraction of excitatory links fexc = 0.5. Insets: Output distribution of the sin- 
gle self-coupled neuron having the target average firing rate /i at the same value 
as the neurons in the network. 

3.2 Noisy autapse 

Let us consider the case of a noisy autapse, when 

x{t) ^ wyit) + ^(t) . (10) 

The neuron receives, beside the autaptic signal ~ wy{t), a random input from 
an external source ~ ^{t) (e.g. from some other neurons in the network). The 
non-autaptic component of the input is drawn from a Gaussian distribution, ^ ~ 
A^(0, cr), where A'^(0, a) denotes a normal distribution with zero mean and variance 

The external input hence perturbs the signal coming from the autapse. From 

10 



0.50 - 



0.25 • 



0.8 

0.-6- 

go.4 

0.2 



- - cj = 



o.oo„ 



• 




k « < 








• 













0.0 



0.1 



0.5 



0.2 
a 



0.6 



fexc 



0.7 



0.3 



0.4 



Figure 4: The Kullback-Leibler (KL) divergence Dx{a,b) of the neuronal firing- 
rate distribution relative to the target exponential distribution (|3]), as a function 
of the standard deviation a of a Gaussian input distribution Px{x), in the presence 
of an autapse (green dots, w = 1 in Eq. ffTOj) ) and in the absence of the autapse 
(red dashed line, w = in Eq. ffTOj) ). Target mean firing rate fi is set to 0.3 . 
Inset: Mean KL divergence (-Da) (see section H]) of the random recurrent neural 
network with = 1000 neurons and mean target firing rate fi = 0.3, as a function 
of the fraction of excitatory links fexc- Note that increase of fexc can be related to 
the decrease in the noisy component of the input that each neuron receives (see 
section 112]) • 

Fig. [2] it is quite obvious that the output distribution of a self-coupled neuron with 
^ = 0, viz. noiseless autapse, in Eq. ffTOj) . deviates substantially from the target 
exponential distribution. The output distribution in the case of noiseless autapse 
is presented in the insets of Fig. [31 However, as we increase the magnitude of 
the external signal, the output distribution of the neuron approaches the opti- 
mal distribution and the KL divergence decreases toward a minimum, see Fig. |H 
Obviously, when a ^ u, the external input dominates, and the two cases with 
and without an autapse become equivalent. Nevertheless, even small amounts of 
noise, that is a < u are sufficient to disrupt the oscillatory behavior of the output 
activity. This happens because of the existence of the second stable fixpoint, for 
certain values of a and b in y = ga,b{y)- As the standard deviation a increases, 
the probability that the firing-rate y will transit toward the second fixpoint also 
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Figure 5: (Left) Mean Kullback-Leibler (KL) divergence (-Da) (color coded) as a 
function of the fraction of excitatory links fexc and target mean firing rate /i in 
a recurrent network of = 1000 neurons and connectivity K = 100. (Right) 
(Dx) as a function of connectivity K and target mean firing rate /i, with excita- 
tory/inhibitory neurons (Top) and projections (Bottom). Fraction of excitatory 
links fexc = 0.8 in both cases. Density plot was evaluated as a linear interpolation 
of the experimentally obtained values represented with green dots. 

increases. Thus, at a certain levels of noise, the activity stochastically escapes in 
short time intervals from the stable fixpoints, and the regular oscillatory behavior 
is destroyed. This also implies that a certain level of decorrelation between the 
input current and output activity has to be reached, if the firing-rate distribution 
is to come as close as possible to the desired target distribution. 

4 Recurrent neural network 

We studied numerically random recurrent neural networks (RRNN) of poly- 
homeostatically adapting neurons where each neuron receives input from K 
pre-synaptic neurons. 

In a first step we consider networks of dual neurons, i.e. a single neuron can 
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have both excitatory and inhibitory projections. In such setup, the synaptic input 
that the ith neuron receives is expressed as 

K 



Xi{t) 



The synaptic weights are selected as Wij = il/y/K, with a probabihty fexc for 
the weight to be positive. The learning rate in ([7]) is set to e = 0.01. We consider 
homogeneous networks where all neurons have identical target average firing rate 
II, determined through (jl]), for the target output distributions Px{y), see Eq. 

For a neuron to ideally map an arbitrary input to the exponential distribution 
with specified mean, it needs a transfer function which can take any functional 
form during the adaption process. This is obviously not the case for the logistic 
function which has only two adaptable parameters. As a result of this limited flex- 
ibility of the transfer function, even when the input of the neuron is independent 



from the output (Fig. Hired dashed line), th e output distribution will neve r idea. 



match the target exponential distribution ( iMarkovic et al. 



2010 



Trieschl . I2OO5I ). 



ly 



4.1 Dynamical behaviors 

To examine the mean deviation of the output activity from the target exponential 
distribution, we have estimated the KL divergence D\{ai,bi), see Eq. ([5]), for 
all neurons in a network of = 1000 neurons, and averaged over the entire 
network and over n = 50 random network realizations. The obtained mean (-Da) is 
presented in Fig. [5] as a function of target mean firing rate fi, fraction of excitatory 
links fexc and network connectivity K. Note that (-Da) is low for high target firing- 
rates and for balanced excitation/inhibition, or for low connectivity K. 

We have observed three distinct dynamical regimes, a pure chaotic regime 
characterized by low values of (-Da), a synchronised oscillatory regime observed 
in random networks with dominating excitatory connections and a intermittent- 
bursting regime observed for balanced excitation/inhibition and small fi. 

To illustrate the difference between these dynamical behaviors, we present in 
Fig. E] the average neural activity (yellow line) and the activity patterns of a ran- 
domly selected neuron in the network of = 1000 units. As we vary the fraction 
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Figure 6: Activity of one, randomly chosen, neuron from the network of N=1000 
neurons with connectivity K=100, depending on the fraction of excitatory hnks 
fexc (top; fi = 0.3) and the mean target firing rate fi (bottom; f^xc = 0.5), where 
the yellow line represents average network activity. The right ordinate shows the 
corresponding mean Kullback-Leibler divergence (see section H]). 

of excitatory links fexc and keep /i fixed (Fig. [6] Top) the dynamics shifts from 
the chaotic phase into the phase of synchronised oscillations. On the other hand, 
reducing the target mean firing rate for balanced excitation/inhibition (Fig. [6] Bot- 
tom) leads to a manifestation of bursts of chaotic activity alternated by periods of 
nearly constant activity. On the right side of the graphs we give the values of the 
corresponding mean Kullback-Leibler divergence (-Da)- In addition, from the oscil- 
latory regime it is possible to transit back into a chaotic or intermittent-bursting 
regime (depending on the value of fi) by reducing the network connectivity K. 
This can be easily seen from the similarity of the density plots on the right and 
left hand side in Fig. [51 
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In Fig. [3] we give an example of the output distributions for two neurons, with 
highest and lowest values of when the dynamics of the neural network is set in 
chaotic regime. The two output distribution are compared with a corresponding 
target exponential distribution. 

Alternatively to randomly selecting a single link as excitatory or inhibitory, 
one can consider a case when a single neuron is selected as either excitatory or 
inhibitory. Thus, all the projections from one neuron are of the same type. We 
have shown (-Da) for such case on the upper right graph of Fig. |5l The absence 
of a visible difference between the upper and the lower graph indicates that the 
intrinsic adaptation, in the low K limit, leads to the same dynamical behavior 
indipendent from having excitatory and inhibitory neurons separeted or neurons 
with both types of projections. 

4.2 Oscillatory behavior 

The network dynamics makes a transition into a synchronized oscillatory regime 
(see Fig. [6]), as /e^c, the fraction of excitatory links, is increased. To better un- 
derstand this oscillatory behavior let us recall the discussion from the previous 
section. In the case of a single self-coupled neuron we showed how a certain 
level of decorrelation, between the output activity and the input signal has to be 
achieved in order for a neuron activity to properly match the target distribution 
(see Fig. H]). The same argument holds in the case of RRNN. Thus, when the 
input is uncorrelated with its output activity (corresponding to the fexc ^ 0.5), 
the output distribution Pa,b{y) closely matches the target distribution px{y). 

The total synaptic input a neuron receives can be divided into two compo- 
nents. The first component is correlated with its own output activity via excita- 
tory recurrent connections. The second component corresponds to the noisy and 
uncorrelated part of the input which results from the competition between inhibi- 
tion and excitation. The first correlated part of the input becomes dominant over 
the noisy second contribution as we increase the fraction of excitatory links fexc- 
The activity therefore starts to follow an oscillatory locked-in trajectory for large 
fractions fexc- 

In the inset of Fig. H] we present the change of mean KL divergence as the 
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number of excitatory connections grows, but for a fixed target mean firing rate. 
We can see that (-Da) increases rapidly once the number of excitatory connections 
starts to dominate. Note that the transition between two phases occurs in a 
shghtly different manner compared to the case of the neuron with a noisy autapse. 
One reason for this difference is that input/output correlations will be amplified 
by additional delayed components of an excitatory feedback a neuron receives. 
This reasoning can be shown to hold by simulating a single neuron with delayed 
coupling autapses driven by the input x(t) — )■ Ylt=o ^kVit ~ t^) + Of) ■ 



4.3 Intermittent bursts of chaotic dynamics 

In the second graph (Fig. [6] Bottom), the dynamics enters an intermediate phase, 
characterized by intermittent bursting of chaotic neural activities, as the target 
mean firing rate is decreased. A closer look into the phase space of intrinsic param- 
eters {ai,bi), of the ith neuron, reveals that the intrinsic parameters approach a 
limiting cycle, similar to the case of a neuron with an autapse. During the regime 
of nearly constant activity Aj(t) ^ the gain steadily increases. Once the gains 
of sufficient number of neurons crosses a certain critical value, the activity of the 
entire network shifts into a chaotic regime. The activity during the chaotic regime 
exceeds the target average activity level /i, thus the gains of all neuron are driven 
back to sub-critical values. Even when reducing the learning rates by several or- 
ders of magnitude this intermittent-bursting behavior persists. Nevertheless, when 
considering constant and supercritical gains for all neurons and allowing only the 
respective biases to adapt, we observe a pure chaotic behavior. This change, which 
arises when reducing the number of degrees of freedoms by considering constant 
gains Oj, is not yet fully understood. A one possible cause could be the use of 
"vanilla" gradient which doesn't take i nto account t he curvature of the manifold 



of probability distributions Pa,b{y) (see lAmaril . 119981 ) . and therefore doesn't point 



into direction of maximal change of Dx{a, b). 
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Figure 7: Nonlinear finite time Lyapunov exponent X{So, t) (see section I4.4p at 
the time step r, for tfie various target average firing rates, with N = 1000, K = 
100, fexc = 0.5 and 6o = 10~^ in all cases. Presented curves are the average of 10^ 
random perturbations. Inset: The limit of dynamical predictability Tp defined as 
a time needed for an error to reach 98% of the saturation level. 

4.4 Sensitivity to external perturbations 

We have evaluated the nonlinear finite time Lyapunov exponent (FTLE), which 
measures the short-term growth rate of initial pert urbations without linearization 



of the time evolution equations ( iDing. fc Lil . 120071 ). 



In practice the FTLE is estimated by considering a small perturbation of the 
trajectory z{t) G along a randomly selected direction 6{0) G M^^, and follow- 
ing the deviation of the perturbated trajectory from the reference orbit, that is 
S{t) = z'{t + r) - z{t + r). The FTLE is the obtained as X{z{t), 6o, r) = i In 
where 6{t) = \\S{t)\\ and 6o = ||5(0)|| ^ 1. The FTLE depends on the starting 
point z{t) of the initial perturbation and on the size of the initial displacement 
^o- The mean FTLE X{6o,t), which is independent from the starting point z{t) 
is evaluated by taking the average of the FTLEs over various points along the 
trajectory, 

1 " 

X{6o,r) = (A(f(t),5o,r))i = - V A(i'(t), 5o, r) . 

n ^-^ 

t=i 

The mean FTLE still depends on the initial displacement 6q. If 6q is chosen to be 
very small one observes initially an exponential growth of the perturbation. For 
this time period the mean FTLE is essentially constant and reduces to the largest 
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Lyapunov exponent. The growth of the perturbation eventually enters a nonlinear 
phase, which is maintained until the deviation from the reference orbit reaches a 
saturation value. Note that FTLE {X{z{t) , 6q, t)) is not necessarily positive for 
all r and z(t), implying that an initial deviation can converge back towards the 
reference trajectory. 

When analysing the changes of the dynamical behavior as we reduced /i, while 
keeping /exc constant at 0.5, we found that when /i is in the range which cor- 
responds to small values of (-Da) (see lower right part of Fig. |5]), the dynamical 
behavior is in a pure chaotic dynamical state with the constant part of the mean 
FTLE X{t) in the range [0.05,0.1] . In this phase the FTLE is positive for all r 
and z{t), thus small initial displacements diverge along every point of the orbit 

As the target mean firing rate is decreased down to ^ = 0.2 we observe a kink 
in the FTLE, with a transition to a second linear time development, see Fig. [71 
The manifestation of the kink corresponds to the occurrence of periods of quasi- 
constant activity as seen in the bottom plot of Fig. El This laminar periods are 
characterised by a negative local Lyapunov exponent, that is small perturbations 
are suppressed during the laminar periods, leading however to a growth of the 
perturbation during the periods of chaotic bursting. 

In this intermittent bursting phase the short term chaotic behavior [t < 200) 
describes the repulsion of two initially close trajectories mainly during the bursting 
regime. The long term behavior (t ^ 200) is also chaotic, as a consequence of 
the intermittent chaotic bursts. The change in the growth of perturbations (see 
Fig. [7]) results from the interplay of the distinct characteristic timescales of the 
intrinsic variables (oj and bi) and the firing rates (vi), with the first being slow 
and the later being fast variables (IBoffetta. Giuliani. Paladin, fc Vulpianil . Il998l ) . 
The occurrence of transiently stable periods of activity also leads to an increase in 
the time Tp, measuring the limit of dynamic predictability, as shown in the inset 
of Fig. [71 The duration of predictability Tp is define d here as the time needed for 
a perturbation to reach 98% of the saturation level (jPing et al.l . 120071 ). 

The dependency of the FTLE on the position of a perturbation along the orbit 
in phase space is presented in Fig. [HI We compared the FTLE, that is X{z{t), 5o, r). 
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Figure 8: Position dependent finite time Lyapunov exponent \{z{t),5Q,T) (see 
section 1131) along the orbit. We considered three cases: (bottom) a chaotic phase 
with /i = 0.25 and f^xc = 0.5, (middle) intermittent-bursting phase with /i = 0.15 
and fexc = 0.5, (top) synchronised oscillatory phase with yU = 0.15 and fexc = 0.7. 
The initial displacement 6q, number of neurons and connectivity K were set to 
10~^, 10^ and 10^, respectively. 

as estimated from six different points along the trajectory z{t), for all three dy- 
namical regimes. The FTLE is then evaluated for 500 consecutive timesteps and 
after the 500th timestep a new perturbation is introduced. In the pure chaotic 
dynamical regime the FTLE is positive for all given initial perturbation points. 
While, in intermittent bursting regime one can also notice negative values of the 
FTLE when the perturbation is initiated within the laminar period. In contrast, 
perturbations starting during the periods of bursts leads to a strictly positive 
FTLEs. In the oscillatory regime the FTLE is negative along the orbit and, simi- 
lar to the single neuron case, the trajectory is unstable during the fast transitions 
from low to high activity states (and vice versa), which results in sharp, positive 
valued, spikes in the FTLE. 

Chaotic dynamics is also observed in the non-adapting limit with e — )■ 0, when- 
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ever the static values of are above the critical value. This is in agreement 
with the results of a large- mean field analysis of an analo gous continuous time 
Hopfield network (ISompolinsky. Crisanti. fc Sommerd . Il988l ). Subcritical static at 
lead, on the other hand, to regular dynamics controlled by point attractors. 



5 Discussion 



Our results show that the introduction of intrinsic plasticity in random recurrent 
neural networks results in ongoing and self-sustained neural activities with non 
trivial dynamical states. For large networks we have observed, depending on the 
specified parameters, three self-organized distinct phases. The network parameters 
include the fraction of excitatory connections, the average connectivity and the 
target average firing rate. These results show that non-synaptic adaptation plays 
an important role in the formation of complex patterns of neural activity. 

An important part of the sensory signals an organism receives result from the 
reactions of the environment to the motor actions taken by the organism itself. The 
complexity of this portion of sensory inputs will then depend on the complexity 
of the organism's own behavior. This consideration indicates that self-generated 
and autonomously sustained neural activity is important for the generation of 
non-trivial behavioral patterns. It is hence more likely that an animal will start 
an explorative behavior if the brain, supporting the body, is able to maintain a 
change in sensory input. In other words, if the dynamics of a neural controller 
would approach, in the absence of sensory inputs, a state of stable and constant 
activity, it would be capable of generating only trivial motor actions. 

As mentioned in the introduction, synaptic plasticity alone drives the dynamics 
of a recurrent network generically toward a f rozen state, independent of the pres- 
ence or absence of sensory input (ISiri et al.l . 120071 1. Synaptic plasticity is thus in 
general, for non-spiking neurons, not sufficient for achieving self sustained activity, 
a likely essential precondition to complex behavioral patterns. 

The relevance of c r itical brain dynamics for both non-linear sensory processing 



investigated intensively. 



(IKinouchi. fc Copellil . |2006|) and for self-sustained neural c o mput ation is being 



Levina. Herrmann, fc Geisell (l2007l 120091) have demon- 
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strated that critical neural activity can be achieved when the depletion of synaptic 
vesicles is included into the dynamics of membrane potential. Under this setup, 
they observe power law scaling of avalanches formed by the activit y of spiking neu - 



rons, a result in agreement with experimental observations (e.g., IChialvd . l2010l ). 
which are supportive of the notion that the brain works in critical regime. How- 
ever, without any form of adaptation to varying sensory stimuli, neurons would 
perform only trivial computations and criticality, or other complex activity pat- 
terns, would generically not arise. Thus, to properly understand the brain dynam- 



ics and cog nitive processes, one must include various forms of plasticity (jTrieschl . 



2005 



20071). 



Here we showed that intrinsic or non-synaptic plasticity will drive a system 
of recurrently interacting neurons, under certain quite general conditions (net- 
work connectivity, ratio of excitation versus inhibition), towards a chaotic phase. 
Synaptic adaption rules, on the other side, are known to generically drive recurrent 
neural networks into a subcritical or frozen state. Our results hence indicate that 
self-organization of neural network dynamics into a critical regime could occur 
whenever intrinsic and synaptic plasticity are both present and relevant. Critical 
neural dynamics would then result from the interplay between synaptic and non 
synaptic adaption processes. 



An analogous line of arguments has been brought forward by 



Per. Hesse. &: Martins 



(120061 ) ■ by demonstrating the relevance of self-organized criticality for the emer- 
gence of exploratory behavior in autonomous agents. Optimal predictability of 
the sens ory-motor cycle is achieved when the neural controller works in a critical 



regime 



Deretal 



20061). We believe that self-organized criticality in biologically 



inspired autonomous recurrent neural networks will exhibit similar patterns of 
complex behavior. Certainly, the complex behavior should persist when including 
the interaction between an agent and the environment, as we plan to do in future 
investigation. 
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